Novel methods for medicinal dosage determination and diagnosis

ABSTRACT

The present invention is directed to systems and methods for diagnosing tissue abnormality or diseases, determining effective drug dosages, and monitoring therapeutic drug treatments. The methods and systems described utilize tissue imaging in situ and computer modeling.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. patent application Ser. No.13/890,198, filed May 8, 2013, and U.S. Pat. No. 8,463,552, filed Oct.18, 2006, and provisional patent application Ser. No. 60/728,034, filedOct. 18, 2005, all of which are hereby incorporated by reference intheir entirety.

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY SPONSOREDRESEARCH AND DEVELOPMENT

This invention was made with government support under grant number RO1NS44598 awarded by the National Institutes of Health (NIH). Thegovernment has certain rights in the invention.

SEQUENCE LISTING OR PROGRAM

A working model of the computer program described below is attached inthe Appendix filed as a computer listing text file.

FIELD OF THE INVENTION

The present invention is directed to systems and methods for diagnosingtissue diseases, determining effective drug dosages, and monitoringtherapeutic drug treatments. The methods and systems described utilizetissue imaging in situ and computer modeling.

BACKGROUND OF THE INVENTION

A large number of psychiatric (i.e. schizophrenia), neurological (i.e.Parkinson's disease), and neurodegenerative (i.e. Huntington's chorea)pathologies involve changes of mental states or conditions based uponchanges in neurotransmitter and receptor balances. Detection of suchchanges may allow for diagnosis well ahead of manifestation of severeclinical symptoms, and knowledge of the nature and the extent of suchchanges is of paramount importance for the determination of therapy. Forinstance, in Parkinson's disease the chronic use of L-DOPA therapy leadsto a progressive diminution in its efficacy. Thus, one would like to beable to monitor the progression of the disease more closely to effectpossible changes in dosing. Similar problems present for many of thecurrently used dopaminergic ligands its schizophrenia. Determination ofthe effects of these therapies upon the brain is very difficult at thepresent time.

The brain provides especially illustrative examples because of variousdifficulties in accessing the tissue directly during life. In a broadersense, though, many biological questions reduce to “what is thesensitivity of this tissue, or this subject, to a given drug or to othermanipulation of a given receptor?”

Some exiting methods allow identification of a region of tissue that isespecially sensitive to a given pharmacological agent (“mapping”), butdo not characterize that response across a range of doses or quantifythe sensitivity. This has uses but usually fails for comparisons acrosslong stretches of time or between subject groups.

Two methodologies have been widely used for the determination of changesin neurotransmitter and receptor dynamics in vivo. These two techniques(Positron Emission Tomography and Single Photon Emission ComputedTomography, PET and SPECT) involve the use of radioactivity. PositronEmission Tomography is a very versatile technique which has been usedsuccessfully for the mapping and sometimes quantification of CerebralBlood Flow (CBF), cerebral glucose metabolism (using[¹⁸F]fluorodexyglucose, FDG) or receptor activity (using radioactivepharmacological ligands), while SPECT has until recently been morelimited to the detection of nonspecific processes. Unfortunately, bothtechniques suffer from severe limitations in spatial and temperalresolution, and cannot generally be proposed for repeated applications.Measuring receptor sensitivity quantitatively with PET or SPECT oftenrequires invasive steps such as arterial catheterization. Moreover, PETis characterized by limited availability and high costs, which arepartly due to the short half-life of many of the cyclotron-producedradiopharmaceuticals which have to be administered.

A number of existing strategies allow quantification of pharmacologicalsensitivity in vivo. Traditionally this is done by measuring abiological signal after repeated experiments, each performed at aspecific drug dosage, across a wide dose range. The resulting data canbe plotted to show a dose-response curve (see FIG. 1). Signals can bequite varied and include clinical information such as blood pressure orsubjective emotional state, endocrine responses, measurement ofchemicals in tissues or body fluids, volume of an organ, cell number orshape, staining and microscopic examination, or measurement ofelectrical signals such as electrocardiogram (EKG) or single-neuronelectrical recordings.

The area of pharmacology that describes and interprets such informationis called pharmacodynamics. Traditionally, data from dose-responsecurves are fit to a mathematical model of drug-receptor interactions.These pharmacodynamic (PD) models allow characterization of thedose-response curve's shape, and position along the dose axis, by asmall number of model parameters. Common models used to describedrug-receptor interactions return parameters such as the B_(max) andK_(D) often used m tracer kinetic studies,, or the E_(max) and EC₅₀ orED₅₀, more commonly reported in clinical studies. The parameter E_(max)refers to the maximum possible magnitude of a measured effect even athigh doses of drug. The concentration of drug in blood or another bodyfluid that produces an effect half the magnitude of E_(max) is denotedEC₅₀. Instead of blood concentration one can characterize theadministered dose of drug similarly; a dose of size ED₅₀ produces aneffect of magnitude E_(max)/2.

All these solutions to calculating quantitative pharmacodynamicinformation are limited in important ways. These include (a) spatiallimitations (such as with needle electrophysiological recordings,microdialysls, or measurements of systemic, e.g. clinical or endocrine,responses; (b) limitations to interpretation at one cellular level, e.g.receptor quantification by PET does not provide a measurement of thesensitivity of the whole complex of receptor and postsynaptic cell(s);and (c) limited pools of subjects, e.g. with autopsy, biopsy, PET,microdialysis, or intrasurgical recordings. Further, PET or SPECTmeasurements of specific receptors can proceed only after extensivetesting of each specific radioligands, including tests of safety, drugdistribution including penetration of the organ of interest,radiodosimetry, displaceability, and pharmacological specificity.

An alternative approach has been called pharmacological challengeimaging. In this approach a receptor-effector system is manipulated byadministration to the organism of a specific pharmaceutical agent, suchas a synthetic agonist at a receptor for an endogenous transmitter orhormone. Additionally, an imaging device records a response of tissue tothis agent. Generally in this case the response, such as glucose oroxygen metabolism, or blood flow, is a general reflection of tissueactivity and can also be observed in response to varied otherphysiological stimuli (such as sensory stimulation, mental work, orphysical exercise). Pharmacological challenge imaging has severalpotential specific advantages: the entire system of receptor andeffector mechanisms is probed simultaneously, responses can be measuredin numerous locations simultaneously (e.g., across the entire brain),and responses can be seen in numerous regions all of which may beaffected by receptor stimulation in the whole organism.

Pharmacologic challenge imaging began with ex vivo rodent smiles withautoradiographic studies of local blood flow or metabolism in tissueslices. In the past several decades PET and SPECT have also been used toobserve changes in humans and other animals. However, only rarely havequantitative measures been used or a wide range of dosses administeredso that one could adequately characterize dose-response curve to thedrug.

More recently, functional magnetic resonance imagine (fMRI) has beenused as the imaging method in pharmacological challenge imaging. Thisapplication has been dubbed pharmacological Magnetic Resonance Imaging(phMRI). The most common fMRI technique is sensitive to changes in bloodoxygen level (Blood Oxygen Level Dependent signal, BOLD). Unfortunatelymost BOLD-sensitive fMRI is nonquantitative. This implies that responsesto drug cannot be determined over a wide range of drug doses.Additionally, quantitative comparisons of drug response over time orbetween groups in an experiment are largely impossible. The situation issomewhat analogous to measuring distance with a stretchable ruler—it mayhe useful over short periods of time or in certain clever experimentaldesigns, but one would view with suspicion an apparent differenceobtained by measurements from the stretchable ruler taken in twodifferent settings.

Other investigators have sought to circumvent the “stretchable ruler”problem of BOLD-sensitive fMRI by adding additional experimental data.Some groups have used what could be called a pure pharmacokinetic (PK)approach. In this approach, blood levels of drug after a single dose ofdrug are computed repeatedly over time from a subject or from a group ofsubjects. The resulting time-concentration curves are used to constraina search for imaging signals that correspond in time so the expectedblood concentrations of the dose of drug. The general plan is to observechanges that occur rapidly in comparison to the to the more gradualartifactual changes in the nonquantitative signal: by analogy, tomeasure quickly, before the ruler stretches too much. One seriousproblem with this method is that it assumes that tissue response is alinear function of drug when in fact two low doses may haveindistinguishable or negligible effect and two higher doses may eachhave maximal effect (FIG. 1). Another problem with this method is thatonly drugs with rapid onset and rapid removal from the system aresuitable. While such drugs exist, e.g. injected nicotine or inhaledcocaine, the realty is that many prescription drugs of interest areselected especially not to have these characteristics. Patients preferto take only one dose of drug a day, so slow drug elimination is oftenpreferred. Since many side effects are concentration-related, slowerabsorption and consequently lower peak concentrations are also desired.

Another attempt to circumvent limitations of BOLD fMRI forpharmacological imaging involves what might be called purepharmacodynamic phMRI. In this method a clinical parameter, such as painrelief by an opiate or a feeling of “high” induced by cocaine, ismeasured repeatedly over time either in the scanner or during a separateexperimental session. The resulting time-response curve is used toconstrain a search for imaging signals that correspond in time to theexpected clinical effect of the drug. Safety or ethical issues aside,this approach suffers from specific problems including (a) the imagingsignal observed cannot easily be interpreted as being a direct effect ofdrug, rather than an indirect effect of pain relief or feeling “high”;(b) only drugs with an easily detectable clinical effect can be used;and (c) studies of a medicament used for a specific disease will usuallyrequire studies only in patients, rather than in control groups, sincethe control subjects have by definition normal function not expected toimprove detectably with administration of drug.

Newer imaging methods may overcome some of the disadvantages ofBOLD-sensitive fMRI for pharmacological challenge imaging. Arterial spinlabeling (ASL) is one such method that promises to provide morequantitative measures. Images taken after administration of MRI signalenhancing agents (“contrast”) can be made quantitative under somecircumstances. These methods have their own disadvantages includingrelatively poor temporal sensitivity or the need to administer thecontrast agent with potential side effects.

However, even if she specific disadvantages of existing imaging methodsfor pharmacological challenge imaging are overcome, they do not generatequantitative pharmacodynamic parameters in current implementations. Infact, in order to generate such quantitative data, all currentpharmacological activation imaging approaches require repeatedexperiments, each after a different dose of a drug of interest, togenerate accurate dose-response curves or otherwise to estimatequantitative pharmacodynamic parameters such as EC50. These repeatedexperiments must also include experiments done at relatively high dosesof the challenge drug (e.g., those that produce at least 90% of themaximal possible signal, or doses that can be many times higher than thedose that produces 50% of the maximal possible signal).

To summarize, the prior art is characterized by a state in whichquantitative measures of pharmacodynamic parameters require one or moreof the following undesirable characteristics; (a) invasive measurement;(b) many measurement (e.g. scanning) sessions; (c) a quantitativemeasure, e.g. a quantitative imaging method; (d) knows pharmacokinetics;(e) observable clinical reactions such as pain relief or a drug “high.”

Without quantitative information some types of clinical information andresearch goals cannot be achieved. The dose of drug needed by a givenpatient cannot be estimated except to say that a response is or is notseen to a given test dose. Diagnostic tests cannot be done unless thedetection threshold for any response to drug happens fortuitously tocorrespond to the optimal drug dose response to which best separates illfrom healthy tissue. Any response to an intervention, other than an “allor none” response, can not be detected. Progression of disease severityover time cannot be monitored.

On the contrary, the present application acknowledges the lack ofexisting methods to provide quantitative pharmacodynamic information onresponses to drugs and simultaneously map those responses in differentregion of an organ or organism. It also acknowledges the lack ofsatisfactory existing quantitative methods for diagnosis or drug dosagedetermination by imaging. The present application uses specificpharmacological stimuli and non-invasive imaging techniques with highspatial resolution and gives a solution to said research and medicalneeds.

SUMMARY OF THE INVENTION

The present invention provides quantitative pharmacodynamic info from asingle imaging session with input images that do not need to bequantitative, in response to modest doses of drug. Certain embodimentsof the present invention may be used to collect information from normalhealthy volunteers, in a very short space of time. Such information istypically useful in designing first (and lengthy) efficacy studies inpatients. One example would be identification in a month, from a dozennormal volunteers with one morning's study each, without psychiatricexpertise, of the likely dose range to test in an initial phase II2-month trial of treatment of human volunteers with major depression.

Certain embodiments of the methods and systems described herein include:

-   -   1) Collect repeated functional imaging data over time (e.g.,        BOLD fMRI) of a body function that in some body part changes in        response to drug (e.g., midbrain changes in blood flow after an        acute dose of levodopa). Perform appropriate motion correction,        artifact reduction, filtering, and (if needed) quantification.    -   2) Give multiple doses of the drug during the time frame imaged,        preferably via a rapidly absorbed route. The doses could be        equal as we have tested in prior work, or could be unequal (e.g.        0.001mg, 0.01mg, 0.1mg).    -   3) As each location imaged, fit a        pharmacokinetic/pharmacodynamic (PK/PD) model to the imaging        time-activity curve over the total scanning period, including as        either knowns or (constrained or unconstrained) variables both        pharmacokinetic information (e.g., distribution halflife,        elimination halflife) and pharmacodynamic information (e.g.,        relative EC₅₀, Hill coefficient). This is described in detail,        below.    -   4) Identify effect locations, organs or volumes of interests        (VOIs) that have the most statistically significant responses,        compared to a null hypothesis such as no change, or change that        fits a polynomial or other function of time that is unrelated to        the experimental design. For example, this may be performed        using a ratio of the summed squared error (SSE) after fitting        the complex PK/PD model to the SSE alter fitting a polynomial        with an equal number of parameters to be fit, the ratio can be        interpreted as an F statistic. Other possible approaches are        obvious to those of ordinary skill in the art.    -   5) Correction of statistical results for multiple comparisons if        multiple volume elements (voxels) arc imaged. For example, this        may be performed using the methods of SPM2, but other methods        are obvious to those of ordinary skill in the art, including        ascertainment of typical statistical distributions in biological        data and judging significance from this actual distribution.        Such methods may better account for the actual distribution in        the statistical images produced in the previous step.    -   6) Extract relative quantitative phrmacodynamic information from        the modeling results at each organ, voxel or VOI.    -   (7) If the administered dose is modeled, the results may be        quantified in terms of ED₅₀ (half-maximally effective dose) or        EC₅₀ relative to the (unknown) concentration of drug in a body        fluid at a given moment relative to administration of the drug.        Optionally, by adding measured concentration(s) of drug in a        body fluid, the absolute EC₅₀ or similar quantitative        pharmacodynamic information may be calculated. Pharmacokinetic        information also may be computed from imaging data or, for        example, from blood concentrations only.    -   (8) The quantitative pharmacodynamic information thus obtained        may then be applied to, e.g., dose determination, diagnosis, or        treatment planning.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a typical example of a dose-response curve from anunrelated study.

FIG. 2 is an overview of the pharmacokinetic-pharmacodynamic (PK-PD)model. It includes a time concentration curve for apomorphine in plasmausing default parameters (top), with arrows indicating the EC₅₀(concentration that produces half the maximal effect).

FIG. 2 (middle) illustrates the shape of the sigmoid Emax model for thedose response curve on a logarithmic scale.

FIG. 2 (bottom, a) is the BOLD v. time scatter diagram used indetermining the EC₅₀ indicated by arrow “a” in FIG. 2 (top).

FIG. 2 (bottom, b) is the BOLD v. time scatter diagram used indetermining the EC₅₀ indicated by arrow “b” in FIG. 2 (top).

FIG. 2 (bottom, c) is the BOLD v. time scatter diagram used indetermining the EC₅₀ indicated by arrow “c” in FIG. 2 (top).

FIG. 3 is an overview of the computer program in the form of a flowchart.

FIG. 4 is a flow chart describing the process of fitting input data to amodel.

FIG. 5 is a flow chart describing the process of finding best-fitparameters for the input data.

FIG. 6 here is a screen shot of the program running, with the setupscreen in the background and the status to the foreground. The programcan be configured to fit each parameter, or iterate through any or allof them as is necessary to practice the invention.

FIG. 7 provides an example of a simulated BOLD time course, both ideal(solid line) and with increasing noise (dots).

FIG. 8 shows predicted sensitivity from validation data, varying bynoise and input EC₅₀.

FIG. 9 also shows predicted sensitivity from validation data, but with amore stringent statistical threshold.

FIG. 10 shows the percentage of voxels within range of the correct inputEC₅₀ for varying noise levels.

FIG. 11 shows the positive predicted value (PPV) for varying noiselevels and input EC₅₀.

FIG. 12 shows a time interactivity curve for one volume of interest inone subject.

FIG. 13 shows three positive predicted value corves for three noiselevels that are similar to the noise levels in the time course of FIG.9.

FIG. 14 shows one slice through an example EC50 image (right), and acorresponding anatomy image (left).

DETAILED DESCRIPTION

I. Introduction

The present invention provides methods and systems that are effective indetermining effective drug dosage, diagnosing diseases associated wishbody organs, including cancers and the monitoring of treatment of suchdiseases. This is accomplished through systems that include applicationof repeated smaller doses of a pharmaceutically active compounds over abrief time interval of minutes such as 2, 4, 5, 10, 15, 20, 30, 45, 60,or 90 minutes, or hours such as 2, 4, 6 8 or 10 hours. The applicationof a plurality of doses is accompanied by computer modeling allowing thepractitioner to calculate relative EC₅₀ for the pharmaceutically activecompound at any point in a target organ. By further sampling a bodilyfluid for drug concentration at any time point including a single timepoint, the practitioner may calculate a EC₅₀ (e.g., 92 ng/mL).

Thus the present invention provides (1) computation of apharmacodynamics parameter (EC₅₀) along with pharmacokinetic parameterspoint in brain, (2) multiple doses of drug, and (3) the computer programthat does the computation. Moreover, certain embodiments of the presentinvention allow collection of information in healthy individuals and isparticularly useful in situations where drug effects may take anextended time to manifest themselves.

II Imaging Devices

Imaging devices suitable for use in the present invention are known inthe art as any imaging device capable of detecting a visible or physicalchange in a bodily organ in situ. Examples include PET, SPECT, andfunctional MRI (fMRI) methods and devices. Visible changes in an organinclude necrosis including scarring, color variations, size variationsand the like. Physical changes include changes in temperature, density,etc. . . .

The present invention is suitable for use in monitoring any organ of abody including, but not limited to, heart, liver, muscle, pancreas,ovaries, testicles, bone, cartilage, eyes, kidneys, and preferablybrain. Identifiable portions of an organ, or effect sites may also beimaged either separately or in conjunction with multiple effect sites,the organ as a whole, or even multiple organs. For purposes of thepresent invention, the term “organ” includes effect sites of the organunless otherwise stated.

III. Acceptable Pharmaceutically-Active Compounds

Compounds suitable for use with the present invention may be anypharmaceutically-active compound that elicits a pharmacodynamicresponse. Such compounds may be formulated into pharmaceuticalcompositions using techniques well known to those in the art. Suitablepharmaceutically-acceptable carriers are available to those in the art;for example, see Remington's Pharmaceutical Sciences, 16^(th) Edition,1980, Mack Publishing Co., (Oslo et al., eds.).

As used herein the language “pharmaceutically acceptable carrier” isintended to include any and all solvents, dispersion media, coatings,antibacterial and antifungal agents, isotonic and absorption delayingagents, and the like, compatible with pharmaceutical administration anddo not interfere with the pharmacokinetic or pharmacodynamic propertiesof the pharmaceutically-active compound. The use of such media andagents for pharmaceutically active substances is well known in the art.Except insofar as any conventional media or agent is incompatible withthe active compound, such media can be used in the compositions of theinvention.

A pharmaceutical composition of the invention is formulated to becompatible with its intended route of administration. Examples of routesof administration include parenteral, e.g., intravenous, intradermal,subcutaneous, oral (e.g., inhalation), transdermal (topical),transmucosal, and rectal administration. Solutions or suspensions usedfor parenteral, intradermal, or subcutaneous application can include thefollowing components: a sterile diluent such as water for injection,saline solution, fixed oils, polyethylene glycols, glycerine, propyleneglycol or other synthetic solvents; antibacterial agents such as benzylalcohol or methyl parabens; antioxidants such as ascorbic acid or sodiumbisulfite; chelating agents such as ethylenediaminetetraacetic acid;buffers such as acetates, citrates or phosphates and agents for theadjustment of tonicity such as sodium chloride or dextrose. pH can beadjusted with acids or bases, such as hydrochloric acid or sodiumhydroxide. The parenteral preparation can be enclosed in ampules,disposable syringes or multiple dose vials made of glass or plastic.

Pharmaceutical compositions suitable for injectable use include sterileaqueous solutions (where water soluble) or dispersions and sterilepowders for the extemporaneous preparation of sterile injectablesolutions or dispersion. For intravenous administration, suitablecarriers include physiological saline, bacteriostatic water, CremophorEL™ (BASF, Parsippany, N.J.) or phosphate buffered saline (PBS).

The carrier can be a solvent or dispersion medium containing, forexample, water, ethanol, polyol (for example, glycerol, propyleneglycol, and liquid polyethylene glycol, and the like), and suitablemixtures thereof. The proper fluidity can be maintained, for example, bythe use of a coating such as lecithin, by the maintenance of therequired particle size in the case of dispersion and by the use ofsurfactants.

Sterile injectable solutions can be prepared by incorporating the activecompound (e.g., a receptor protein or anti-receptor antibody) in therequired amount in an appropriate solvent with one or a combination ofingredients enumerated above, as required, followed by filteredsterilization. Generally, dispersions are prepared by incorporating theactive compound into a sterile vehicle that contains a basic dispersionmedium and the required other ingredients from those enumerated above.In the case of sterile powders for the preparation of sterile injectablesolutions, the preferred methods of preparation are vacuum drying andfreeze-drying which yields a powder of the active ingredient plus anyadditional desired ingredient from a previously sterile-filteredsolution thereof.

Oral compositions generally include an inert diluent or an ediblecarrier. They can be enclosed in gelatin capsules or into tablets. Fororal administration, the agent can be contained in enteric forms tosurvive the stomach or further coated or mixed to be released in aparticular region of the GI tract by known methods. For the purpose oforal therapeutic administration, the active compound can be incorporatedwith excipients and used in the form of tablets, troches, or capsules.Oral compositions can also be prepared using a fluid carrier for use asa mouthwash, wherein the compound in the fluid carrier is applied orallyand swished and expectorated or swallowed. Pharmaceutically compatiblebinding agents, and/or adjuvant materials can be included as part of thecomposition. The tablets, pills, capsules, troches and the like cancontain any of the following ingredients, or compounds of a similarnature: a binder such as microcrystalline cellulose, gum tragacanth orgelatin; an excipient such as starch or lactose, a disintegrating agentsuch as alginic acid, Primogel, or corn starch; a lubricant such asmagnesium stearate or Sterotes; a glidant such as colloidal silicondioxide; a sweetening agent such as sucrose or saccharin; or a flavoringagent such as peppermint, methyl salicylate, or orange flavoring.

For administration by inhalation, the compounds are delivered in theform of an aerosol spray from pressured container or dispenser whichcontains a suitable propellant, e.g., a gas such as carbon dioxide, or anebulizer.

Systemic administration can also be by transmucosal or transdermalmeans. For transmucosal or transdermal administration, penetrantsappropriate to the barrier to be permeated are used in the formulation.Such penetrants are generally known in the art, and include, forexample, for transmucosal administration, detergents, bile salts, andfusidic acid derivatives. Transmucosal administration can beaccomplished through the use of nasal sprays or suppositories. Fortransdermal administration, the active compounds are formulated intoointments, salves, gels, or creams as generally known in the art.

The compounds can also be prepared in the form of suppositories (e.g.,with conventional suppository bases such as cocoa butter and otherglycerides) or retention enemas for rectal delivery.

It is especially advantageous to formulate oral or parenteralcomposition in dosage unit form for ease of administration anduniformity of dosage. “Dosage unit form” as used herein refers tophysically discrete units suited as unitary dosages for the subject;each unit containing a predetermined quantity of active compoundcalculated to produce the desired therapeutic effect in association withthe required pharmaceutical carrier. The specification for the dosageunit forms of the invention are dictated by and directly dependent onthe unique characteristics of the active compound and the particularpharmacological effect to be achieved, and the limitations inherent inthe art of compounding such an active compound for the treatment ofindividuals.

IV. PK-PD Model and Statistical Methods

Theory. The novel approach depends in part on the recognition that themost accurately measured variable in most imaging experiments is time.By giving repeated doses of drug and measuring responses to each doseover time intervals short enough to minimize time-dependent artifactualsignal drift, one can compute a quantitative measure of sensitivity todrug in a single imaging session, even with a nonquantitative imagingmethod. This process is summarized in FIG. 2. In the upper right panel,a typical time concentration curve for a drug is plotted usingtraditional pharmacokinetic modeling (solid curve). The three dashedlines represent three different subjects (or three different tissueregions) with varying sensitivity to drug. Line a represents a regionwith high sensitivity (low EC₅₀), whereas lines b and e representregions with moderate and low sensitivity. To estimate the expectedtime-response functions of the various subjects (or regions) requiresthe composition of the time-concentration curve with theconcentration-response curve shown in the middle right graph in FIG. 2.As shown in the lower right panels of FIG. 2, the resultant tissuetime-response curves are markedly different. For sensitive subject orregion a, the first dose of drug produces a maximal response. Forrelatively insensitive region e, only the later doses have a substantialeffect. The following paragraphs formalize this concept.

The pharmacokinetic-pharmacodynamic (PK-PD) model describes therelationship between the dose of a specific drug and the correspondingeffect. This is diagrammatically depicted in FIG. 2 (See Holford,Nicholas H. G., Sheiner, Lewis B. (1982) Kinetics of PharmacologicResponse. Pharmacol Ther. 16: 143-166). The pharmacokinetic aspect ofthe model relates the dose of the drug to the concentration of the drugin a compartment, e.g., in plasma or some effect site. Thepharmacodynamic aspect of the model relates the concentration of drug ina compartment to the drug effect. The combination of these aspects isthe PK/PD model (see FIG. 1).

The invention may be practiced with any number of doses being applied toan individual and/or effect site. At its logical extreme, a single longadministration of drug is equivalent to numerous small, brief doses ofdrug given without interval. By way of an example, four essentiallyinstantaneous IV infusions may be given during a 40-minute phMRI BOLDscan, spaced eight minutes apart beginning 8 minutes after the start ofthe scanning session. The PK model for these infusions can be describedby a simple one-compartment model or by more complex models discussed inany standard pharmacology text. A sigmoid E_(max) model may be used todescribe the PD model. The model is characterized by the maximal effectof drug at high doses, E_(max), the Hill coefficient or sigmoidicityindex n, which models the number of drug molecules required to activatethe receptor and can be understood as describing the steepness of thesigmoid curve at its inflection point, and EC₅₀, which for the case n=1is the drug concentration that produces an effect half as large asE_(max). The middle right graph in FIG. 2 shows this curve on alogarithmic x axis. The input to the curve is plasms concentration of adrug, and the effect that is measured and modeled is the imaging signal.

In the equations below, z(t) describes the plasma concentration of drugwhen K doses of drug, D_(n), are given at times t_(n), and u(t) is theunit step function. The model also includes a fixed time delay (shift)t_(x) and a halflife t_(half) for loss of drug from plasma. This plasmaconcentration z(t) then becomes the input to the sigmoidconcentration-effect model. To add realism the model also adds drift ofsignal unrelated to drug input, modeled by a polynomial of degree M,poly_(M)(t). The sum comprises the model for the imaging signal, herecalled tissue_(model)(t).

The equations that describe the model are given below:

${{poly}_{M}(t)} = {\sum\limits_{n = 0}^{M}{a_{n}*t^{n}}}$${z(t)} = {\sum\limits_{n = 1}^{K}{D_{n}*0.5^{(\frac{t - t_{s} - t_{n}}{t - {half}})}*{u\left( {t - t_{s} - t_{n}} \right)}}}$${{tissue}_{model}(t)} = {\frac{E_{\max}*\left( {z(t)} \right)^{n}}{\left( {EC}_{50} \right)^{n} + \left( {z(t)} \right)^{n}} + {{poly}_{M}(t)}}$tissue_(noise)(t) = poly_(N)(t)

To best fit the model to the data, we originally used theLevenberg-Marquardt method. This method basically tries to minimize thedifference between the real data for each effect site and the givenmodel function by using two other methods for curve fitting, first usingthe steepest descent method when far away from the minimum error, andthe inverse-Hessian method when getting closer to the minimum error(See: Press, William H., Teukolsky, Saul A., Vetterling, William T.,Flannery, Brian P. (1988) Numerical Recipes in C++, 2^(nd) Ed.,Cambridge University Press, Cambridge, UK, pg 688.). The currentimplementation of the model fitting uses a combination of iteration overselected allowed values in parameter space and ordinary least squaresfitting by matrix inversion and multiplication, a method described insection VII below. Other possible parameter estimation methods areobvious to those with ordinary skill in the art.

V. Determining Dosages

The present invention provides methods for determining the effectivenessof pharmaceutical active compounds. For example, the present inventionprovides rapid methods for determining an effective dose for apharmaceutically active compound, and is capable of being practiced inthis regard on healthy individuals. In another aspect, the presentinvention allows rapid patient screening to determine whether aparticular drug or drug regime (cocktail of two or more pharmaceuticallyactive compounds) will be effective in treating individual patients.E.g., the present invention may be used to determine whichanticonvulsant would be effective and at what estimated for a givenpatient. Finding correct effective medications and dosages areespecially critical where the medication is expensive, toxic, oressential for preserving life. Rapid determination of correct effectivedosage is also important where under- or overshooting the dosage wouldbe dangerous or where clinical effect would take weeks to monitor (e.g.antidepressants).

In the present invention, effective medications and dosages aredetermined by administering a plurality of doses of pharmaceuticallyactive compound to a patient. The patient may be a healthy individual ora sick one. Time course imaging of a target organ or effect site is thenperformed and alterations in a delectable property of the target organmonitored. From these observations an EC₅₀ for the administeredpharmaceutically active compound may be determined as described hereinand an effective dosage calculated.

A second exemplary approach involves tracking therapeutic progressduring treatment. For example, Parkinson's disease treatment usuallyinvolves administration of L-DOPA, the precursor of the neurotransmitterDopamine (DA), sometimes combined with substances that block peripheralconversion of L-DOPA to DA. Although L-DOPA replacement therapy hasproven to be very efficacious, it is difficult to optimize the dosage.Too high dosing may lead in the short term to dyskinesia, while in thelong term it may even accelerate the degradation and loss of DA neurons.Optimization of the dose by neurological examination alone is atime-consuming and possibly dangerous approach. The present inventionallows monitoring of dosage and effect in a manner that circumvents suchdifficulties. This is described in detail in the example below. Briefly,through short term monitoring of brain voxel sites known to respond toL-DOPA, an effective dose for the drug may be determined without riskingdangerous over or under dosing or delaying effective treatment.Individual screening is an important aspect of this embodiment of theinvention as it is certain that individual patients will show adifferent response to the various drugs, and it may be impossible orimpractical to determine the optimal treatment for each patient byclinical examination alone. The present invention offers a quantitative,objective end efficient way to compare the efficacy of various drugs inthe treatment of PD.

An additional example is provided by the problem of choosing anappropriate first dose(s) for an early Phase II/III treatment study ofmajor depression. Although tolerability can be estimated from Phase Istudies, estimating the approximate population efficacious dose frompreclinical studies or from normal humans is fraught with difficulty. Aclinical study requires at least 6-8 weeks per subject and often takesover a year to complete. Estimation in healthy or depressed subjects, ina single imaging session, of the EC₅₀ of a metabolic or hemodynamicregional brain response to modest (below EC₅₀) doses of drug by themethod described above could potentially shorten the development cyclesubstantially. Furthermore, if healthy and depressed subjects havesimilar acute brain tissue imaging responses to drug, determining theapproximate dose for an initial Phase II study could be determinedwithout the need to recruit ill subjects and without psychiatricexpertise for subject selection and response monitoring.

Again such dosage determinations are not limited to any particulartissue. For example, a vasodilator drug may be assessed by imagingcoronary artery diameter in vivo during drug infusion.

VI. Diagnosis

The present invention may also be used in detection, diagnosis andguidance of therapy involving a diseased organ, and is particularlyuseful in neurological, neurodegenerative, and psychiatric diseasediagnosis. Such diseases include but are not limited to, Parkinson'sdisease (PD), Alzheimer's disease (AD), Huntington's chorea, andschizophrenia. By way of example, PD is associated with loss ofdopaminergic neurons and is known to remain asymptomatic for a prolongedperiod. Often clinical symptoms are only observed after 70% of thedopaminergic innervation has been lost. By measuring voxel changes, forexample, in she brain in response to a dopaminergic agonist such asamphetamine, asymptomatic individuals may be identified at a stage ofdisease progression where treatment may substantially extend thepatient's quality of life. The feasibility of such an approach isdemonstrated in the example below in human subjects.

As another illustrative example, the method can equally well be appliedto diagnosis, staging and guidance of therapy of Alzheimer's diseasewhich affects four million people in the US alone. As already mentioned,successful application of PET is hampered by the high cost and therequirements for cyclotron and technical support. Measurements of rCBFby SPECT, on the other hand, have been found to be “of limitedapplication for identifying mildly demented AD patients, . . . and oflimited value for distinguishing moderately to severely dementedpatients with AD or vascular dementia” (27). Although success oftherapeutic treatments of AD is, at present, still limited, earlydiagnosis is nevertheless important (a) in order to rule out othercauses like tumors or stroke which might be treated successfully and (b)because the few therapeutic agents known today work best in the earlystages of the disease.

When supplied with the disclosure provided herein, it will be obvious toone of skill in the art that the applications to diagnosis describedherein may be applied in an analogous manner to any tissue where adrug/tissue reaction combination is known in the art and associated witha disease state of the tissue responsive to the drug.

VII. Computer Method

The computer method determines quantitative pharmacodynamic parametersfrom biological data when such data are acquired during apharmacological activation experiment involving administration ofvarious doses of drug (see FIG. 3 and 4 for overview). As built, thismethod can be applied to any time-series biological data, though it isoptimized for data from an imaging system. The method involves thecombination of several components described below. For brevity, andbecause they are not essential characteristics of the method, obvioussteps such as tiding data from a disk file into memory or writingresults to disk are omitted.

The program includes a module that predicts the blood concentration overtime, C(t), in response to a supplied description of the times andquantities of administered doses of drug and a standard one-compartmentpharmacokinetic model as described in section IV above (see FIG. 5). Thepredicted function C(t) also includes a variable delay to account fornonzero administration times or noninstantaneous absorption, or topartially account for hysteresis.

The program also predicts the response of the tissue or organism tovarying blood concentrations of drug, E(C), using a typical sigmoidpharmacodynamic model as described in section IV above.

The composition of these two functions, E(C(t)), varies depending uponthe supplied dose administration details and the values supplied for thevarious pharmacokinetic and pharmacodynamic model parameters describedabove, such as Emax, EC₅₀, the Hill coefficient n, the eliminationhalf-life, the time delay, and so on.

The program also allows the user to select which model parameters willbe held constant and at what values, and what will be the startingvalues and allowable range of the remaining (fit) parameters.

To compare two possible one-dimensional functions of the time variable,as sampled discretely over the time interval given, the program computesan error function, either the chi square error function or the summedsquared error. This error or cost function is used to measure thegoodness of fit of the prediction E(C(t)) to the actual suppliedbiological data, with better goodness of fit being reflected by lowererror.

The program samples parameter space and minimizes error by adjusting theparameter values within their allowable domains using any of threemethods. An early version of the program used a nonlinear curve fittingalgorithm as described in Section IV above. Subsequent applications ofthis program to simulation data revealed that this method did notperform better than the method described next. For all parameters butEmax and the polynomial terms (“iterated parameters”), the programiterates over a discrete set of values representative of eachparameter's allowed domain. Approximately eight to 10 values provideadequate sampling for each parameter but this could of course be varied.Using these parameters the predicted shape of E(C(t)) is computed usingE_(max)=1. For convenience call this result E₁(C(t)), The program thendetermines the value of E_(max) and of the polynomial coefficients a,that optimally fit the predicted full PK-PD model

a ₀+a ₁ t+. . . +a _(n) t ^(n)+E_(max)·E₁(C(t))

to the provided biological data. This is done exactly by the method ofordinary least squares using matrix inversion and multiplication;algorithms for these methods are given in many standard numericalmethods tests. The cost function is then computed for each set ofparameter values. The set of values for the iterated parameters thatminimized the error function, over all allowed combinations of theseparameters, is identified as the optimal set of values and returned tothe user. Other methods of optimizing the iterated parameters will beobvious to those with ordinary skill in the art.

The program also computes a statistic that quantifies the goodness offit of the full PK-PD model to the provided biological data. Severalpossible methods for statistical model comparison have been described.The current implementation uses the following approach. The number ofPK-PD model parameters chosen to be fit to the data by the user is addedto the degree of the polynomial described in the preceding paragraph.Call the sum m. A polynomial of degree m is optimally fit to thesupplied data, again using least squares methods to provide the exactsolution. Note that this polynomial has the same number of freeparameters (m) as the full PK-PD model. The squared error at each timepoint represented in the provided data, summed over all such timepoints, is summed for the full PK-PD model computed as described in thepreceding paragraph and also for the best-fit polynomial of degree m.The ratio (summed squared error for polynomial)/(summed squared errortor best-fit full model) is reported to the user and is interpreted asan F statistic with m and [(number of time points sampled in theprovided biological data)−m−2] degrees of freedom.

The computer program supplied in the appendix (on CD-ROM) is a workingmodel of this method (see FIG. 6). It was written and compiled usingMicrosoft Visual Basic and is operable on a computer using a Pentium orcompatible processor and the Microsoft Windows XP operating system with.NET version 2.0. The program prompts for image files supplied in the“.4dfp” format, which is a widely used floating point ANALYZE-compatiblefour-dimensional image “multivolume” with details represented in theInterfile Header Format (.IFH) files, ANALYZE header files, and imagesavailable at http://www.purl.org/net/kbmd/b2k/. Other options set by theuser are obvious to those with ordinary skill in the art.

VIII. Systems

The present invention also provides systems for carrying out the methodsdescribed herein. Typically such systems include an imaging device forobserving the tissue under study and a computing device for performingthe calculations necessary to model the system and/or determineeffective dosages for pharmaceutically active compounds under study. Theappendix to this application provides exemplary computer code forcarrying out such calculations, for example on a personal computer. Thepresent invention contemplates other implementations of the methoddescribed herein, and as such the computer code presented in theappendix should be considered exemplary only and not limiting on thepresent invention, for example, other computer languages and platformsmay be employed. In addition, other arrangements of procedures andfunctions may be used to accomplish the present invention.

EXAMPLE 1 Introduction

We tested the performance of the methods and system described in thisapplication using simulated data.

Predicted time-imaging curves tissue_(model)(t) were generated for awide range of plausible values for EC₅₀, half-life, and other PK-PDmodel parameters. A value of E_(max)=10 was used for all simulations.Gaussian noise was added to tissue_(model)(t) to better reflectreal-world situations. A random noise function of time was generated1000 times for each level of noise (quantified as the standard deviationof the noise for comparison to the known E_(max)) to allow testing ofsensitivity and specificity of the curve-fitting methods.

FIG. 7 shows examples of tissue_(model)(t) (solid line) and the sum oftissue_(model)(t) and noise (dots) for different levels of noise. It isapparent that at high levels of noise, curve fitting will be difficult.

To test how likely the curve-fitting procedure was to report a fit todata in the absence of an (intentional) signal, noise was alsorepeatedly generated and added to a polynomial of low degree. Withoutnoise the algorithm could not fit another curve better than thepolynomial and F should always be low. With substantial noise, however,the resulting data may be fit better by another curve such as thatgenerated by the model.

Parameter estimation

Parameters were estimated using the computer method described above inthis specification. Data first were filtered temporally to minimizenoise, replacing the data at each time point with the median of the dataacquired over a longer interval centered on that time point. For thedata described here a 45-second interval was used.

For the simulation testing described below, we divided predeterminedranges of each PK-PD model parameter into approximately 10 levels eachto simplify testing. In other words, the program could produce onlycertain answers (e.g., EC₅₀ could be only a set of specified valuesbetween 0.1 and 10). These values were appropriate to the range of inputdata supplied the program during simulation testing. This quantizationwas felt to be reasonable since in practice precision is most importantif the EC₅₀ is within an order of magnitude or so of the achieved bloodconcentration; values outside that range tend to have similar minimaleffect or similar maximal effect. For testing, the desired (“correct”)answer was chosen to be a value not available to the program, e.g. 0.43,so as not to favorably bias the results. For the data shown here,halflife was fixed (not fit) and the Hill coefficient n was set to 1.

Statistical test of goodness of model fit to stimulated data

Summed squared error across all time points at a given voxel was used toquantify how well the data at the voxel, voxel(t), were fit by a giventime-activity curve tissue_(model)(t) generated by the PK-PD model and agiven set of model parameters. A comparison function that did notincorporate any information about drug administration was used tostatistically test how well the model fit the data. For simplicity thecomparison function chosen was a polynomial with the same number ofdegrees of freedom as the PK-PD model tissue_(model)(t). Specifically,if k PK-PD model parameters were used to generate z(t), to which apolynomial of degree M, poly_(m)(t), was added to givetissue_(model)(t), the comparison function was poly_(N)(t), where N=M+k.The ratio F of the summed squared error for poly_(N)(t) to the summedsquared error for tissue_(model)(t) is computed and saved to create astatistical image reflecting the improvement in the fit to the data byincorporating knowledge of drug administration times and the PK-PDmodel. Higher values for F indicate better fit for the model.Probability that the model fit better than the comparison polynomial bychance was computed by interpreting F as an F test statistic with k and([number of time points in voxel(t)]−k−2) degrees of freedom.

Simulation testing: test statistics

For each combination of parameters tested, the model was fit to 1000voxels at each noise level. The following summary statistics wereobtained for each set of parameters and noise level:

-   -   mean and range for F and each of the fit model parameters (EC₅₀,        n, etc.)    -   sensitivity=fraction of times (of 1000) that the primary        parameter of interest (EC₅₀) returned by the program was        “correct” i.e. if the value of EC₅₀ used to generate the data        was 0.43, and the possible answers included {0.10, 0.50, 1.0, .        . . }, then only the nearest possible answers, 0.10 and 0.50,        were counted as correct. In other words, how often does the        program return as close as possible to the desired value?    -   positive predictive value (PPV) was defined as follows. A given        allowed output value of EC₅₀ could have been returned by the        program from input generated from a different EC₅₀ value pins        noise. The PPV was computed with denominator equal to the number        of times that the given EC₅₀ value was output by the program        across all sets of 1000 voxels generated by all sets of input        model parameters, and numerator equal to the number of those        times when the output was the “correct” answer for the input. In        other words, for an EC₅₀ output value of 0.5, of all the times        the program output EC₅₀=0.5, what fraction came from an input        EC₅₀ value nearer 0.5 than any other output value?

It should be noted that PPV depends on the prior probability, i.e. howoften the specified output value was supposed to be produced, so theactual PPV in a given experimental situation may differ from thatcomputed here. However, the PPV can be computed from prior probability,sensitivity and specificity when these are known. The first two of theseare described above. Specificity was computed from the(polynomial+noise) data, and was taken to be 1−p, where p is thelikelihood (fraction of times in 1000 attempts) that the programreturned a given value for EC₅₀. Specificity for F was defined similarlyas one minus the likelihood in the (polynomial+noise) data that Fexceeded a given threshold.

Results: model fitting

FIG. 8 shows how often the model fit the curve significantly better thandid the polynomial comparison function. As visible in this figure, themodel fits better essentially all the time unless the noise level isquite high. Even in the presence of substantial noise, the model alsofits better for sensitive regions (i.e. low EC₅₀).

FIG. 9 uses a more stringent statistical threshold for F correspondingto the Bonferroni correction for a typical number of voxels in a brainimaging experiment (64,000). The same descriptive results obtain.

Specificity: The model never produced a significantly better fit to thenull hypothesis (line plus noise) data, i.e. specificity=1.0 for thisdata set.

Results: corrections of EC₅₀ estimate

FIG. 10 shows on the vertical axis the percent of voxels in range of thedesired (input) value for EC₅₀, as a function of input EC₅₀ and noiselevel. The program found the correct value 100% of the time for lowlevels of noise or for low EC₅₀ values. The sensitivity remained highexcept at high levels of noise or for voxels with low sensitivity.

FIG. 11 shows positive predictive value (PPV) on the vertical axis. Thisfigure demonstrates that, of all input time courses for which thecomputer program returned a given value for EC₅₀, most were in factgenerated using that value for EC₅₀. Positive predictive values were100% for regions with low noise or high sensitivity. A very similarpattern was seen when examining only voxels for which the model fitsubstantially better than the comparison polynomial (data not shown).

Summary

If the modest assumptions of the method are met, the novel methodperforms very well even with modest levels of noise. Confidence isjustifiably high that regions identified as being of high sensitivity ofdrug are in fact so. The assumptions are reasonable and are iteratedhere. (a) The sigmoid pharmacodynamic model is appropriate for a givendrug and imaging measures. (b) Artifactual signal drift is reasonablyapproximated by a low-degree polynomial. (c) Measurement errorapproximates a normal curve.

EXAMPLE 2

To determine how the noise levels in the simulation data just describedcorrespond to realistic biological data, we acquired data usingBOLD-sensitive fMRI imaging of an anesthetized monkey withadministration of the dopamine agonist SKF82958.

Experiments were approved in advance by the Washington University AnimalStudies Committee and due care was exercised to prevent discomfort tothe animal. During imaging sessions the animal's health was monitored byveterinary staff from the Washington University Division of ComparativeMedicine. Two male M. fascicularis monkeys were studied.

Imaging was done while the animals were ventilated with 1.5% isofluranein oxygen.

We obtained repeated whole-brain BOLD-sensitive fMRI images at arepetition interval of one every 3 seconds over a 40-minute period.

At 10 minutes into the imaging session a steady-rate, slow intravenousinfusion of the dopamine D1 agonist SKF82958 was commenced finishing at15 minutes, with a total administered dose of 0.1 mg/kg. We havepreviously shown that this dose of drug is safe and does not affectwhole-brain mean quantitative blood flow, hot produces reliable effectson regional blood flow consistent with local neuronal activation (Blacket al. J Neurophysiol 2000; 84(1):549-557).

Volumes of Interest (VOIs) were traced on structural brain imagesacquired at the same session. These included left and right caudate,left end right putamen, striatum collectively, whole brain, and amidbrain region. Time-activity curves from each VOI were extracted andsubmitted to the computer method described above.

A time-activity curve for one VOI in one subject is shown in FIG. 12.From this figure it is seen that E_(max) with this experimental systemis of magnitude at least 20.

A straight line was fit by the method of least squares to the VOI datafrom the first 15 minutes of the experiment, i.e. before any drug wasadministered. The standard deviation (S.D.) of the residual signal aftersubtracting this line from the data was computed. Across different VOIsand the two subjects, this standard deviation varied from 0.21 to 3.33.The S.D. for all striatal regions was in the range 0.97-2.09.

Since E_(max)≥20, the ratio S.D./E_(max)≤S.D./20. This ratio wascomputed and compared to S.D./E_(max)=S.D./10 for the simulationexperiment. Thus in the worst case (E_(max)=20), a typical striatalregion with S.D.˜1.5 corresponds to the fourth lowest noise level in thesimulation data, which had S.D./E_(max)=0.75/10. The positive predictivevalue for this level of noise, depending on the input EC₅₀ value, isshown as the middle of the three curves on FIG. 13. The upper and lowercurves on this figure correspond to S.D./E_(max)=0.5/10 andS.D./E_(max)=1.0/10 respectively.

In summary, using real biological data, the noise level and effect sizeseen with an appropriate pharmacological stimulus, combined with theresults of the simulation testing data, suggest that there is highconfidence in the EC₅₀ values returned for regions identified as havinghigh sensitivity to drug (low EC₅₀) by the methods and systems describedin this application.

EXAMPLE 3 Introduction

This example describes the implementation and testing of apharmacokinetic-pharmacodynamic model for analysts of functional MRIdata. The data illustrates mapping and quantifying a brain BOLD responseto an acute dose of apomorphine in living humans.

Many important neurological and psychiatric illnesses involveabnormalities of dopamine function. An in vivo method so measuredopaminergic system sensitivity would benefit diagnosis and research.Positron emission tomography (PET) has been used to measure regionalcerebral sensitivity to dopamimetic challenges, but PET is not ideal forlongitudinal studies or for research in children.

Regulatory Approval

Approved by FDA (IND ∩62999, 63000) and WUSM Human Subjects Committee(IRB). Subjects gave written informed consent.

Study Population

6 normal control and 9 Parkinson's disease (PD) patients were selectedto participate in the study.

Study Day Overview

The PD patients attended one MRI day on which they were given anapomorphine challenge. The normal control subjects took part in two MRIvisits: an apomorphine challenge and a saline control day. Domperidonewas given at 20 mg p.o. q.i.d. for 6 doses, with the last dose given30-60 minutes prior to the administration of apomorphine. IV catheterswere placed in arm veins for infusion of apomorphine (total dose, 10μg/kg) or a saline placebo. Apomorphine was given during a 40-min fMRIsession, as described below. Afterward, subjects wore asked to rate anyside effects on visual analog scales.

MR Methods

Anatomic images and BOLD sequences were obtained using a 1.5 T SiemensVisions system equipped with boosted EPI gradients. Asymmetric spin echoprovided BOLD-sensitive images [Hershey et al. Biol Psychiatry 2004;55(9):916-925].

Pharmacologic Challenge

Four i.v. infusions were given during the 40-minute phMRI BOLD scan.Each consisted of 0.0025 mg/kg apomorphine given over 30 seconds. Theywere given at 8, 16, 24, and 32 minutes into the fMRI session (see topof FIG. 2),

PK-PD Model

Apomorphine pharmacokinetics were described by a 2-compartment linearmodel with variable parameters. FIG. 2 (top) shows a time/concentrationcurve for apomorphine in plasma using default parameters. We used asigmoid concentration-response curve to describe pharmacodynamics(middle of FIG. 2). FIG. 2 (bottom) shows the predicted BOLD response toapomorphine, plus gaussian noise, given three different values for EC₅₀(concentration that produces half maximal effect).

The following parameters were estimated from individual voxel time: BOLDcurves by an iterative least-squares fit method.

TABLE 1 PK/PD parameters fit using custom software Parameter Startingvalue Ratio of EC₅₀ to first plasma peak 2.0 Elimination halflife ofapomorphine 41.0 min Time delay from plasma peaks to BOLD peaks 0 secLine of polynomial coefficents, 0 a₁ (in a₀ + a₁x + a₂x² + a₃x³) Peakdifference between data and polynomial (E_(max)) 0

Image Analysis

We examined two models of brain BOLD signal under the null hypothesis ofno drug-related activation; (1) straight line drift; (2) cubicpolynomial. Each was compared with the best-fit curve described fey thevariables in Table 1. The ratio F of summed squared error was used toquantify variance better explained by the full model than by lite lineor polynomial.

Results

Analysis of Apomorphine Scans

A typical EC₅₀ image is shown in FIG. 14. Testing continues at this dateto determine how accurately this method differentiates healthy subjectsfrom Parkinson disease patterns and how consistently high-sensitivity(low EC₅₀) regions are identified across subjects.

Conclusions

We have been able to implement a PD/PK model for analysis of fMRI data.By giving repeated small doses of drug, the timing of any brain responsecan be expected to depend on the sensitivity of that brain region todrug. For instance, a region with high sensitivity to drug may show amaximal response beginning with the first dose, whereas a region withlower sensitivity to drug would show a response only after repeateddoses, or not at all. In other words, this method potentially convertsan unknown pharmacologic parameter (EC₅₀) to a measurable variable (timecourse of BOLD signal).

Since the blood arriving at all parts of brain has essentially identicalconcentration of the test drug, the relative sensitivity of differentregions of brain can be directly compared. If additionally one canmeasure the concentration of the drug in blood, the relativesensitivities can be quantified, e.g. the EC₅₀ (concentration thatproduces half-maximal effect) can be computed and given in ng/mL. Themethod could be applied to many other drug challenges. Futurerefinements may include higher doses of apomorphine, bettercurve-fitting, or more quantitative MR methods such as arterial spinlabeling.

ADDITIONAL EMBODIMENTS

The methodology of any of the methods or the system discussed above maybe further enhanced in various ways. One such way would be to addexplicit hysteresis modeling to the dose response of thepharmaceutically active compound by adding an effect site delaycompartment to the pharmacokinetic-pharmacodynamic model. Such anenhancement would take into account the dose response effects of drugssuch as levodopa, whose biological effect at a given moment reflects notonly the current blood level but also recent blood levels. Additionalimprovements whose implementation will be straightforward may includemore complex pharmacokinetic modeling, e.g. a two-compartment model withdistribution and elimination modeled separately, or more complexpharmacodynamic modeling. Long single-dose infusions can be modeled bythese methods and system and are distinct from prior art in the methodused for analysts. Additionally, equal timing and amount of doses is notrequired by the methods or system taught here. The computer method isnot limited in application to imaging data or drug-effect modeling butcould equally characterize drug-receptor binding (K_(D)) or perhapsenzyme kinetics.

Although the foregoing invention has been described in some detail byway of illustration and example for clarity and understanding, it willbe readily apparent to one of ordinary skill in the art in light of theteachings of this invention that certain changes and modifications maybe made thereto without departing from the spirit and scope of theappended claims.

All publications and patent applications cited in this specification areherein incorporated by reference as if each individual publication orpatent application were specifically and individually indicated to beincorporated by reference.

APPENDIX

The APPENDIX filed herein as a Computer Listing text file contains thefollowing files in IBM-PC format compatible with the MS-Windowsoperating system, created Jan. 5th, 2007, containing the followingMicrosoft Visual Basic .NET source code files (© 2006 Kevin Black andJonathan Koller):

File Name File size (bytes) Date of creation ApoFit.Designervb 642711/29/2005 14:21 ApoFit.vb 55314 11/29/2005 14:21 DesignDose,Designervb7804 10/2/2006 9:55 DesignDose.vb 164 10/2/2006 9:55 FitSettings. vb1220 11/29/2005 14:21 FramesToSkipDialog.Designer.vb 2801 4/19/200613:00 FramesToSkipDialog.vb 1428 4/19/2006 13:00 Image4dfp.vb 1389111/29/2005 14:21 LogFile.vb 485 8/7/2006 9:06 MainWindow.Designervb61463 11/29/2005 14:21 MainWindow.vb 8171 11/29/2005 14:21 MatOps.vb3923 11/30/2005 11:27 MLFit,vb 17933 11/29/2005 14:21 ParamSettings.vb1487 11/29/2005 14:21 Stats.vb 17578 11/29/2005 14:21

1. A method of planning treatment for a patient with Parkinson'sDisease, the method comprising: a) administering to the patient withParkinson's disease, one or more doses of levodopa; b) collectingbiological imaging data of one or more brain regions of the patient withParkinson's disease in response to the one or more doses of levodopaadministered in a), over a period of time where the response to levodopahas not reached a steady state, and c) repeatedly measuring levodopaconcentrations in a body fluid from the patient with Parkinson'sdisease, over the period of time, d) fitting predeterminedpharmacokinetic-pharmacodynamic models to the information from step b)and step c), to determine constant ke, or t½ eq; e) comparing theconstant ke, or t½ eq, to historical values from healthy humans, and, orpatients with Parkinson's disease, and f) wherein subsequent treatmentof the patient with Parkinson's disease is determined based on thecomparison in e).
 2. The method of claim 1, wherein the historicalvalues are drawn from said patient with Parkinson's disease.
 3. Themethod of claim 1, wherein the imaging data is collected on a magneticresonance imaging device.
 4. The method of claim 1, wherein the imagingdata is collected on a magnetic resonance imaging device sensitive to ablood oxygen level dependent signal.
 5. The method of claim 1, whereinthe imaging data is collected using Positron Emission Tomography.
 6. Themethod of claim 1, wherein the imaging data is collected on a SinglePhoton Emission Computed Tomography imaging device.
 7. The method ofclaim 1, wherein the imaging data consists of perfusion imagingcollected on a magnetic resonance imaging device.
 8. A method forimproving a measurement of quantitative pharmacodynamic parameters forlevodopa therapy in a human with Parkinson's disease, the methodcomprising: a) administering to the human with Parkinson's disease asingle dose of levodopa; b) collecting biological imaging data of thehuman with Parkinson's disease's response to the single dose of levodopaadministered in a), the biological imaging data comprising, brainactivity images before, during, and after the single dose of levodopa isadministered, c) measuring levodopa plasma concentrations, d) fittingpredetermined pharmacokinetic-pharmacodynamic models to said datacollected in b) and c) to determine the quantitative pharmacodynamicparameters for said single dose of levodopa, comprising: i) setting upfit parameters and choosing a fit method; ii) computing best-fitparameters and their goodness of fit to the biological data using theselected fit method; iii) computing the statistical significance of thebest-fit parameters; and iv) generating an output. v) comparing theoutput to historical output from healthy humans, and, or patients withParkinson's disease.
 9. The method of claim 8, wherein the historicaloutput is drawn from said patient with Parkinson's disease.
 10. Themethod of claim 8, wherein the imaging data is collected on a magneticresonance imaging device. 11 .The method of claim 8, wherein the imagingdata is collected on a magnetic resonance imaging device sensitive to ablood oxygen level dependent signal.
 12. A method for improving ameasurement of quantitative pharmacodynamic parameters for levodopa in ahuman with Parkinson's disease, the method comprising: a) administeringto the human a single dose of levodopa; b) collecting biological imagingdata of the human with Parkinson's disease response to the single doseof levodopa administered in a), wherein, the biological imaging dataconsists of the response of the human with Parkinson's disease collectedbefore, during and after levodopa is administered and c) repeatedlymeasuring levodopa plasma concentrations d) fitting said predeterminedpharmacokinetic-pharmacodynamic models to said data, collected in b) andc) consisting of calculating an effect compartment rate constant valueKe that collapses the hysteresis curve to a dose-effect function, forthe affected brain regions, and e) comparing the effect compartment rateconstant value Ke that collapses the hysteresis curve to a dose-effectfunction, for the affected brain regions, to the effect compartment rateconstant value Ke that collapses the hysteresis curve to a dose-effectfunction for the affected brain regions in a group of humans sufferingfrom Parkinson's, to determine the relative severity for the individualhuman suffering from Parkinson's disease.
 13. A method of improving ameasurement of Parkinson's disease severity in an individual human, themethod comprising: a) administering to the human with Parkinson'sdisease one or more doses of levodopa; b) repeatedly collectingbiological imaging data from one or more brain regions of the human withParkinson's disease over a period of time during which the response tolevodopa in the one or more brain regions has not reached a steadystate; c) repeatedly measuring levodopa concentrations in a body fluidfrom the patient with Parkinson's disease, over the period of time, d)fitting predetermined pharmacokinetic-pharmacodynamic models to theinformation from step b) and step c), to determine the constants ke ort½ eq, and e) comparing the constant ke or t½ eq to historical values ofke or t½ eq from healthy humans, and, or patients with Parkinson'sdisease.
 14. The method of claim 13, wherein the method is performedwithout data from a clinically measured response.
 15. The method ofclaim 13, wherein the historical values of ke or t½ eq from healthypatients, and or, patients with Parkinson's disease patients, consistsof historical values of ke or t½ eq from said patient with Parkinson'sdisease.
 16. A method of improving a measurement of Parkinson's diseaseseverity in one or more regions in the brain in an individual human withParkinson's disease, the method comprising: a) administering to thehuman with Parkinson's disease one or more doses of levodopa; b)repeatedly collecting biological imaging data of one or more brainregions of the human with Parkinson's disease over a period of timeduring which the response to levodopa has not reached a steady state; c)repeatedly measuring levodopa concentrations in a body fluid from thepatient with Parkinson's disease, over the period of time; d) fittingpredetermined pharmacokinetic-pharmacodynamic models to the informationfrom step b) and step c), to determine the constants ke or t½ eq; and e)comparing the constant ke or t½ eq to known values of ke or t½ eq fromhealthy humans, and, or patients with Parkinson's disease.
 17. Themethod of claim 16, wherein the known values of ke or t½ eq are drawnfrom said patient with Parkinson's disease.
 18. The method of claim 16in which severity of Parkinson's disease is determined simultaneouslyfor more than one brain regions.
 19. A method of improving adetermination of the rate of progression of Parkinson's disease in ahuman with Parkinson's disease, the method comprising: a) performingnon-concurrently, on each of two or more study days: i) administering tothe human with Parkinson's disease one or more doses of levodopa; ii)repeatedly collecting biological imaging data one or more brain regionsof the human with Parkinson's disease over a period of time during whichthe response to levodopa has not reached a steady state; iii) repeatedlymeasuring levodopa concentrations in a body fluid from the patient withParkinson's disease, over the period of time, and b) fittingpredetermined pharmacokinetic-pharmacodynamic models to the informationfrom each study day to determine the constants ke or t½ eq applying tothat study day, and c) computing the rate of progression of Parkinson'sdisease in the said human as the slope of the line that best fits theset of data points from all study days.
 20. A method of improving adetermination of whether a medication or biological therapy slows a rateof Parkinson's disease progression, comprising: a) measuring the rate ofprogression of Parkinson's disease according to the method of claim 39in a group of people with Parkinson's disease who are given themedication or biological therapy between the first and last study days,and b) measuring the rate of progression of Parkinson's diseaseaccording to the method of claim 39 in a group of people withParkinson's disease who are not given the medication or biologicaltherapy between the first and last study days, and c) comparing therates of progression of Parkinson's disease between the two groups usingconventional statistical tests.